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We analyse systematically, from the viewpoint of the nonlinear physics of solitary waves, the 
effect of the spatial dimension {D — 1,2,3) on the structure and stability of the Bose-Einstein 
condensates (BECs) trapped in an external anisotropic parabolic potential. While for the positive 
scattering length the stationary ground-state solutions of the Gross-Pitaevskii equation are shown 
' to be always stable independently of the spatial dimension, for the negative scattering length the 

^\ , ground-state condensate is stable only in the ID and 2D cases, whereas in the 3D case it becomes 

On ' unstable. A direct link between nonlinear modes of BECs and (bright and dark) solitary waves of 

the nonlinear Schrodinger equation is demonstrated for all the dimensions. 



in 
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I. INTRODUCTION 



Experimental realisations of Bose-Einstein condensation (BEC) in ultracold and dilute atomic gases [y^] generated 
extensive theoretical studies of the macroscopic dynamics of condensed atomic clouds (see, e.g., Ref. 3 and references 
^ . therein). In experiments, the BEC atoms are trapped in a three-dimensional, generally anisotropic, external potential 
OO ' created by a magnetic trap, and their collective dynamics can be described by the well-known Gross-Pitaevskii (GP) 
, equation, 

o 

O' in^ = -^S/l^ + V{r)^ + Uo\^\'^, (1) 



dt 2m 



where ^'(r, i) is the macroscopic wave function of a condensate in the _D-dimensional space, V{'r) is a parabolic 
trapping potential, and the parameter Uq = ^Trti'a/m characterises the two-particle interaction proportional to the 
s-wave scaterring length a. When a > 0, the interaction between the particles in the condensate is repulsive, as in 
^ , most current experiments jl], whereas for a < the interaction is attractive 
'■^ ' From the viewpoint of the nonlinear dynamics of solitary waves, the case of the negative scattering length is the 
Ch , most interesting one. Indeed, it is well known that the solutions of Eq. (^ without a parabolic potential display 
^ collapse for both D = 2 and _D = 3, so that no stable spatially localised solutions exist in those dimensions at all. 
A trapping potential allows a stabilisation of otherwise unstable and collapsing solutions and, as has been already 
shown for the 3D case with the help of different approximate techniques |^,^, there exists a critical number N^i- of 
the total condensate particles N, defined through the normalisation of the condensate wavefunction as 



N = J (2) 

such that the condensate is stable for N < N^, and it becomes unstable for N > Ncr, and then its local density j^E'p 
grows to infinity. This critical value of the particles' number is readily observed in experiments 

In this paper we analyse systematically, in the framework of the GP equation (^, the structure and stability of a 
trapped condensate depending on its dimension D. In spite of the fact that the first experimental results were obtained 
for the BEC in a 3D anisotropic trap, the cases of lower dimensions are also of great importance and, moreover, they 
have a clear physical motivation. First of all, Eq. (^ with D = 1 appears as the fundamental model of BEC in highly 
anisotropic traps of the axial symmetry 



V{f) = ^{rl+Xz'), (3) 

provided A = w^/'-^l <C 1, i.e. for cigar-shaped traps. In this case, the transverse structure of the condensate is 
controlled by a trapping potential, and it can be described by the eigenmodes of a two-dimensional isotropic harmonic 
oscillator. Then, after averaging over the cross-section, Eq. (|l|) appears as a ID model with rcnormalised parameters, 
and it describes the longitudinal dynamics of the condensate (see Ref. 6 and Sec. 2 below). 

Secondly, in the case D = 2, the model (||) can be employed to describe the condensate dynamics in the so-called 
pancake traps (A ^ 1), and it can be related to recent experiments with quasi- two-dimensional condensates or coherent 
atomic systems 
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The paper is organised as follows. In Sec. 2 we consider the ID case and find numerically the ground-state solutions 
for the trapped condensate in both a < and a > cases. In the former case, the problem is found to be remarkably 
similar to that of the soliton propagation in optical fibers in the presence of the dispersion management ||^ . The most 
interesting case of a 2D trap is considered in Sec. 3, where we discuss the necessary conditions for the instability of the 
radially symmetric condensate in a trap, and also analyse the invariants of the stationary solutions found numerically. 
We point out, for the first time to our knowledge, that 2D condensates can be stable, in spite of the possibility of 
collapse in the model described by the 2D GP equation. Similar results are briefly presented in Sec. 4 for the 3D 
BECs of radial symmetry, and they provide a rigorous background for the variational results obtained earlier by 
different methods [Q. In all the cases, we also find the condensate structure for the positive scattering length where 
the instability is absent in all the dimensions. 



II. QUASI-ONE-DIMENSIONAL TRAP 
A. Model 

To derive the ID GP equation from Eq. we assume that the parabolic potential V{'r) describes a cigar-shaped 
trap. For the potential (^), this means that A <C 1, and the transverse structure of the condensate, being close to 
a Gaussian in shape, is mostly defined by the trapping potential. We are interested in the stationary solutions of 
the form ^'(r, i) = ^'(r)e^Tr*, where fi has the meaning of the chemical potential. Measuring the spatial variables 
in the units of the longitudinal harmonic oscillator length aho — {h/mujy/xy^^ , and the wavefunction, in units of 
{huj /2Uoy/xy ^'^ , we obtain the following stationary dimensionless equation: 



$ + cr|*P$ = 0, (4) 



where /i' — 2/i-\/A/?iw, (^,?y,C) = {x,y,z)/aho, and the sign a = sgn(a) = ±1 in front of the nonlinear term is defined 
by the sign of the s-wave scattering length of two-body interaction. 

We assume that in Eq. (|lj) the nonlinear interaction is weak relative to the trapping potential force, i.e. Uo/muj ~ 
A ^ 1. Then, it follows from Eq. (^) that the transverse structure of the condensate is of order of A, and the 
condensate has a cigar-like shape. Therefore, we can look for solutions of Eq. (^) in the form, 

$(r,C) ^^WVXC), (5) 



where r = -^/^^ + 77^, and is a solution of the auxiliary problem for the 2D quantum harmonic oscillator 

V±</- + 7(A- A-ij'^^ = 0, (6) 
which we take in the form of the no-node ground state, 

0oM=C^exp(-^), 7= A, (7) 

C being constant. To preserve all the information about the 3D condensate in an asymmetric trap describing its 
properties by the ID GP equation for the longitudinal profile, we impose the normalisation. 



/oo poo 
(j){r)dxdy = 2tt (j){r)rdr 
-00 Jo 



= 1 (8) 



that yields = I/ttn/A. 

After substituting the solution (H) with the function (0) into Eq. (Q), dividing by and integrating over the 
transverse cross-section , 77) of the cigar-shaped condensate, we finally obtain the following ID stationary GP equation 
(see also the similar results in Ref. 6) 



+ /3^_^2^ + ^|^|2^^0, (9) 
where f3 = fi' — j, and the normalisation (H) for 0(r) is used. 
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Importantly, the number of the condensate particles is now defined as A'' = {huj/2Uo^/X)Q, where 

/oo 
l^l'dC (10) 
-oo 

is the integral of motion for the normahsed nonstationary GP equation. 

B. Stationary Ground-State Solutions 

Equation (||) includes all the terms of the same order, and it describes the longitudinal profile of the condensate 
state in an anisotropic trap. In the linear limit, i.e. when formally ct ^ 0, Eq. (^ becomes the well-known equation 
for a harmonic quantum oscillator. Its localised solutions exist only for (3 = /3„ = 1 + 2n (n = 0, 1, 2, . . .), and they 
are defined through the Gauss- Hermite polynomials, ^n{C) = c„e~^ ^^Hn{C), where 

H.iO - i-ireC'^^'-^^. (11) 

so that Hq = 1, Hi = 2C, etc. 

To describe the effect of weak nonlinearity, we can employ a perturbation theory based on the expansion of the 
general solution of Eq. ^ in the infinite set of eigenfunctions defined by Eq. (^ij) . The similar approach has been 
used earlier in the theory of the dispersion- managed optical solitons (see, e.g., Ref. 10), and it allows us to describe 
the effect of weak interaction on the condensate structure and stability. These results will be presented elsewhere . 

In the opposite limit, i.e. when the nonlinear term is much larger than the parabolic potential, localised solutions 
exist only for /? < and a = +1, and they are described by the stationary nonlinear Schrodinger (NLS) equation. 
The one-soliton solution is 



MO = iTi^y (12) 
cosh(Cv— P) 

so that Qs — 4:^/—f3 coincides with the soliton invariant. 

In general, the ground-state solution of Eq. ^ can be found only numerically. Figures 1(a) and 1(b) present 
some examples of the numerically found solutions of Eq. (^) for several values of the dimensionless parameter /5, for 
both negative (a) and positive (b) scattering length. For (3^1, i.e. in the limit of the harmonic oscillator mode, 
the solution is close to Gaussian for both the cases. When (3 deviates from 1, the solution profile is defined by the 
type of nonlinearity. For attraction (a = +1), the profile approaches the sech-type soliton ([T^), whereas for repulsion 
(cr = —1) the solution flattens, cf. Fig. 1(a) and Fig. 1(b). 

In Fig. 2 we show the dependence of the invariant Q on the parameter /3, for both the families of localised solutions, 
corresponding to two different signs of the scattering length. The dotted line shows the limit of the soliton solution of 
the NLS equation without a trapping potential. In the asymptotic, i.e. say for (3 < —2, the curve Qsif3) coincides with 
the invariant Q for the BEG condensate in a trap. This means that for such localised condensates the effect of the 
trap is negligible and the condensate function becomes localised mostly due to an attractive interparticle interaction. 

To discuss the stability of localised solutions, we refer to a number of the corresponding problems well analysed 
in nonlinear optics of guided modes and solitary waves [Q. In terms of the so-called Vakhitov-Kolokolov stability 
criterion, the condition dQ/d{—(3) > gives the solution stability for the case cr = -1-1. Therefore, all ID solitons for 
the attraction case are stable. A detailed discussion of the soliton stability for this case can be found in Ref. 9. 

Formally, the case a = —1 corresponds to the opposite sign of the derivative dQ/d(3. However, this case of the 
so-called defocusing nonlinearity is also well-known in nonlinear optics of guided waves, and it has been analysed for 
different types of the trapping potential. In particular, it is known that the Vakhitov-Kolokolov criterion does not 
apply to this case, and all the family of spatially localised solutions are stable ||l^. Extension of all those results to 
the case of a parabolic trap is a straightforward technical problem. 



C. Higher-Order Modes and Multi-Soliton States 

It is well known that in the linear limit cr ^ 0, Eq. (^ possesses a discrete set of localised modes corresponding to 
the different orders of the Gauss-Hermite polynomials. We demonstrate that all those modes can be readily obtained 
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for the nonhnear problem as well, giving a continuation of the Gauss-Hermite modes to a nonlinear regime. Figure 
3 shows several examples of those modes for both negative {a = +1) and positive {a = —1) scattering length, 
respectively. In the limit /3 — > 1, those modes are defined by the eigenfunctions of the linear harmonic oscillator. The 
effect of nonlincarity is different for the negative and positive scattering length. For the negative scattering length 
(attraction), the higher-order modes transform into multi-soliton states consisting of a sequence of solitary waves with 
alternating phases [see Figs. 3(a) and 3(b)]. This is further confirmed by the analysis of the invariant Q vs. /3, where 
all the branches of the higher-order modes tend asymptotically to the soliton dependences Qn ^ {N + l)Qs, where 
N is the order of the mode {N = 0, 1, . . .). Examples of the dependences Q{/3) for the first- and second-order modes 
are shown in Fig. 4. From the physical point of view, the higher-order modes in the case of attractive interaction 
correspond to a balance between repulsion of out-of-phase bright solitons and attraction of the trapping potential. It 
is clear that such a balance can be easily destroyed by, for example, making the soliton amplitudes different. However, 
the analysis of the stability of such higher-order multi-soliton modes is beyond the scope of this paper. 

For the positive scattering length (a = —1), the higher-order modes transform into a sequence of dark solitons [ p^ , 
so that the first-order mode corresponds to a single dark soliton, the second-order mode to a pair of dark solitons, 
etc. [see Figs. 3(c) and 3(d)]. Again, these stationary solutions satisfy the force balance condition - repulsion between 
dark solitons is exactly compensated by an attractive force of the trapping potential. 



III. QUASI-TWO-DIMENSIONAL TRAP 
A. Model and Stationary Ground-State Solutions 

The case of the two-dimensional (2D) GP equation can be associated with a quasi-condensate of atomic hydrogen 
[0 or quasi-2D gas of laser cooled atoms Derivation of the GP equation from the first principles of scattering 
theory in the two-dimensional geometry is not trivial, and it has been recently shown ||l5|] that the correct form of 
the 2D GP equation (|l|) should have the 2D interaction potential Uq in the form Uo = — (27r?i^/m) In^^(fca), where 
< ka < 1, a is the scattering length, and the characteristic wavenumber k can be approximated as fc 1/aho- Here, 
we derive the 2D GP equation directly from its 3D form of Eq. (|^), for the case of a pancake trap when A 1. 

In the case of a pancake trap, in the parabolic potential (|^) we assume the parameter A is large, i.e. A 3> 1. Then, 
the longitudinal profile of the condensate is controlled by the parabolic potential ~ (mAa;^/2)z^. Measuring the 
spatial variables in the units of the transverse harmonic oscillator length (Ti/tocj)^/^ and the wavefunction in units of 
(?itj/2C/o)^^^, from Eq. (|]) we obtain the following stationary equation, 

+ m'* - [(^^ + V^) + AC^] $ + ^ 0, (13) 

where /x' = 2/x/?ia;. For A 3> 1, the longitudinal structure of the condensate is squeezed in the z-direction and its size 
in this direction is of the order of A^^. Therefore, we can look for stationary solutions in the form, 

$(r,C)=V(0'/'(C), (14) 

where ipi''^) depends on the radial coordinate r = \/ + 2/^i a-nd this time the function (f>{C) is a solution of the ID 
quantum harmonic oscillator. 



'2, 



^^^+7</'-AC'0 = O, (15) 



with the Gaussian form, 



0o(C) = Cexp|^-^^, 7-VA. (16) 

Normalisation condition 0Q(C)dC = 1 yields C — tt^/^'A^^^. 

Substituting the solution°(^^ into Eq. ( p^ ) and averaging over C, we obtain the following 2D equation for ij{r), 

|^ + ;|:) V^ + Z^^-r^^ + dlV-lV^O, (17) 



4 



where a — ±1 and /3 = /i' — 7. Equation (|l7| ) is the 2D GP equation for the stationary solutions without any angular 
dependence. In the linear limit, i.e. when we neglect the nonlinear term (cr ^ 0), the solution exists only at P — 2, 
and it can be written in the form V'o('') = Co cxp (— where Co is defined by normalisation. In the opposite 
limit, the localised solution exists only for a — +1, and it is described by the radially symmetric solitary wave of the 
2D NLS equation. 

We integrate Eq. (|l7| ) numerically and find a family of radially symmetric localised solutions ip{r) for the condensate 
ground state, for both a ~ +1 (attraction) and cr = —1 (repulsion). Some results are presented in Figs. 5(a,b), for 
different values of the parameter /3. The structure of the localised solutions is similar to those in the ID case. However, 
in the limit of large negative /?, the solution transforms into the 2D NLS soliton, known to be a self-similar radially 
symmetric solution corresponding to a critical collapse. This is further illustrated in Fig. 6, where we present the 
dependence of the invariant Q on the parameter /3. The 2D NLS soliton exists only for /3 < and cr = +1, and it 
corresponds to the fixed value of the invariant, Qs ~ 11.7, shown as a dotted straight line in Fig. 6. 



B. Stability and collapse 



Stability of the 2D condensates in a parabolic trap is an important issue. In particular, as was shown by Berge [l^ ] 
and Tsurumi and Wadati the presence of a trapping potential does not remove collapse from the 2D GP equation. 
This result is very much expected because it has been known for a long time that, in the case of the 2D NLS equation 
with a parabolic potential, there exists an exact analytical transformation that allows the potential to be removed 
from the 2D NLS equation |l^. This means that for some classes of the input conditions of the GP equation, the 
collapse should be observed even in the presence of a trapping potential. Nevertheless, as follows from our results 
summarised in Fig. 6, for the attractive case (cr — +1) the derivative dQ/d{—(3) is positive and it does not change 
the sign. Moreover, according to the collapse conditions derived in Ref. 15, for the radius R of the time-dependent 
2D condensate the following equation holds 



where H is the system Hamiltonian 



H ^2n 



and R is defined as 





d^/j ^ 




dr 




= 2tt 



8H - A{R^ 



^dr. 



(18) 



(19) 



(20) 



surely 
tends 



Equation (^8|) means that the condition _ff < is a sufficient condition for collapse, because the value (i?^) 
becomes negative for some value of t and, therefore, the wavefunction becomes singular and collapses as {R^ 
to zero. 

When H is zero or positive, the dynamics depend on the applied perturbation, so that for a large enough perturbation 
we can lower the stationary value of the Hamiltonian H . In Fig. 7 we present the dependences H{(3) and H{Q) for 
the family of stationary solutions found numerically for the attraction case cr = -1-1. It is clear that the whole family 
of stationary solutions corresponds to _ff > 0. Therefore, the ground-state is expected to be stable everywhere, 
but its dynamics should be sensitive to the perturbation amplitude for larger (3 since a perturbation can make the 
Hamiltonian negative. This statement should be further verified by numerical simulations of the 2D GP equation. 



C. Higher-Order Modes and Vortices 

Similar to the ID case discussed above, the higher-order modes of a harmonic oscillator can be found in the 2D 
case too, thus providing an analytical continuation of the 2D Gauss-Hermite modes. The set of such solutions is 
broader because it includes a continuation of a superposition of different modes. For example, the mode tpooir) is the 
fundamental ground-state solution, that continues the zero-order Gauss-Hermite mode ~ Hq{x)Hq{ii), as discussed 
above. In general, the modes 4'nmix,y) have no radial symmetry and they include both dipole-like and vortex- like 
states. 
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The simplest higher-order radially symmetric mode describes the condensate with a single vortex, shown in Fig. 8. 
Figures 8(a) and 8(b) show the density profile of the condensate with a vortex state, for both negative (cr — +1) and 
positive (a = —1) scattering length, respectively. Figure 8(c) presents the dependence of the invariant Q vs. /3 for the 
vortex state in comparison with the ground-state mode. Because the vortex state is a nonlinear mode that extends 
the corresponding mode of the linear system, it approaches the value = 4. 



IV. THREE-DIMENSIONAL TRAP OF RADIAL SYMMETRY 



The case of a radially symmetric 3D trap is the most studied in the literature, so that the corresponding results are 
well known, and therefore we reproduce some of them just for the completeness of the general picture of the stationary 
states. Similar results for the stationary states can be found, e.g. in Refs. 3 and 17, whereas the BEC stability has 
been recently discussed by Huepe et al. . 

Figures 9(a) and 9(b) show the density profiles of the 3D condensate of radial symmetry for different values of j3, 
for both negative (cr = -1-1) and positive {a = —1) scattering length. While in the case of repulsion (cr = —1), the 
condensate broadens for larger (3 remaining always stable, for the attractive interparticle interaction (cr = +1), the 
condensate is stable only for f3 < where /3cr = 0.72 corresponds to the maximum value Qcr ~ 14.45 [see Fig. 
10]. This result is consistent with the Vakhitov-Kolokolov stability criterion dQ/d{—(3) > 0. The critical value /3cr 
corresponds to a critical value of the particles N„ and it has been already observed in experiment Without the 
trapping potential i.e. in the limit of the 3D NLS equation, all localised solutions for cr = +1 are unstable (shown as 
a dotted curve in Fig. 10). 



V. CONCLUSIONS 



We have presented a systematic study of the ground-state and higher-order spatially localised solutions of the GP 
equation for the Bose-Einstein condensates in a parabolic trap of different dimensions. While many results for the 
radially symmetric ground-state and vortex modes in 3D are available in the literature for the condensates with the 
positive scattering length (i.e., for repulsive interaction of the condensate atoms), little is known for the condensates 
with the negative scattering length, especially for the 2D traps. In particular, we have presented the results indicating 
that, in spite of the collapse condition in the 2D case, the family of 2D localised solutions of the GP equation can 
be stable. Additionally, we have clarified the meaning of higher-order localised modes, i.e. nonlinear Gauss-Hermite 
eigenmodes, and demonstrated their link to multi-soliton states. We believe the systematic study of the stationary 
states in all the dimensions carried out here from the viewpoint of the nonlinear physics of localised states and 
solitary waves allows us to close a gap in the literature and deepens our knowledge about the condensate properties 
and stability. 
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a = -1 



FIG. 2. Dependence of the invariant Q, defined by Eq. (jlOfK on the parameter /3 for the family of the ground-state localised 
solutions of Eq. ^ . Dotted - result for the ID NLS soliton p2|) without a trapping potential. 
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FIG. 4. Invariant Q vs. /3 for the first- and second- order nonlinear modes, shown together with the corresponding dependence 
for the fundamental mode of Fig. 2. Dotted lines are given by the values Qs, 2Qs, and SQs, respectively. 





FIG. 6. Dependence of the invariant Q on the parameter (3 for the family of the ground-state solutions of Eq. (y_7|). Dotted 
- the limit of the 2D NLS soliton without a trapping potential. 
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FIG. 7. Dependence of the Hamiltonian H on (a) the parameter j3, and (b) on the invariant Q for the family of 2D radiaUy 
symmetric solutions of Eq. (|l7|) at a = +1. 
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FIG. 8. (a,b) Vortex profiles and (c) the corresponding invariant Q for the ground-state and vortex localised modes of the 
2D GP equation. 
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FIG. 10. Invariant Q vs. the parameter /3 for the family of the radially symmetric ground-state solutions of the 3D GP 
equation. 
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